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1 Introduction 

While QCD factorization methods are well established in the case of scattering observablcs 
involving a single large mass scale pQ, the treatment of processes with multiple mass scales is 
rather more subtle. In multi-scale processes, generalized factorization formulas are required [2] if 
one is to control perturbative large logarithms to higher orders of perturbation theory and to de- 
scribe appropriately nonperturbative physics in the initial and final states of the collision. Such 
generalized factorization formulas typically involve transverse- momentum dependent (TMD), 
or "unintegrated" , parton distribution and parton decay functions [2j [3] . 

A full treatment of factorization at unintegrated level, valid uniformly over the whole space, 
is yet to be achieved [4]. Results however exist which apply in specific phase space regions. One 
such case is given by QCD in the high-energy, or small-x, limit [5], in which factorization of 
TMD gluon distributions holds [BJ, in correspondence to pcrturbativcly-rcsummed coefficient 
functions. Results based on small- x TMD distributions have for instance been used in analyses 
of inclusive observablcs such as deeply inelastic structure functions (see [7] for a review) and in 
Monte Carlo simulations of the exclusive structure of final states in hadronic collisions (see [8]). 

Many aspects of the experimental program at the Large Hadron Collider (LHC) depend on 
the analysis of processes containing multiple hard scales, and will be influenced by improved 
formulations of factorization in QCD at unintegrated level. In this article we give a concise 
overview of recent progress on TMD distributions and their use for simulations of final states in 
hadronic collisions by shower Monte Carlo event generators. Most of TMD computational tools 
have so far been developed within a quenched approximation in which only gluon and valence 
quark effects are taken into account at TMD level. We describe results of recent work [Qj [10] 
to go beyond this approximation by including sea quark effects, and we present first numerical 
applications to forward Z-boson production. 

The article is organized as follows. In Sec. 2 we recall operator definitions for parton 
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distribution functions (PDFs) and basic issues associated with their TMD generalization. In 
Sec. 3 we consider parton branching methods and the role of unintcgrated distributions in shower 
Monte Carlo generators. In Sec. 4 we discuss the approach [9] to incorporate effects from quark 
emission and flavor-singlet sea-quark contributions in the framework of transverse-momentum 
dependent parton showers. We summarize in Sec. 5. 



2 TMD parton distribution functions 

In this section we introduce the parton correlation functions used to define parton distributions. 
In Subsec. 2.1 we consider operator matrix elements for TMD distributions, and discuss current 
open issues, including lightcone divergences and factorization breaking effects. In Subsec. 3.2 
we focus on the case of small x, related factorization results, and introduce applications which 
are the subject of the sections that follow. 



2.1 Operator matrix elements 

The relevance of consistent operator definitions for parton kj_ distributions was emphasized 
long ago in the context of Sudakov processes [TT], jet physics [T2"| . exclusive production }13) . 
spin physics [T3] ■ The approach commonly used to ensure gauge invariance is to generalize the 
coordinate-space matrix elements that define ordinary parton distribution functions (pdfs) |15j 
to the case of field operators at non-lightconc distances. For instance, for the quark distribution 
one has (Fig. [1]) 

f(y) = (P\iJ(y)Vj(n) 7 + V (nW(0)\P) . (1) 

Here ip are the quark fields evaluated at distance y = (0, y~ , y±_), where y± is in general nonzero, 
and V are eikonal-line operators in direction n, 

Vy(n)=Vexp(ig s dr n"A M (y + r n)j , (2) 

which we require to make the matrix element gauge- invariant. The unintcgrated, or TMD, 
quark distribution is obtained from the double Fourier transform in y~ and of /. 




Figure 1: Correlator of two quark fields at distance y. 

While Eq. ([1]) works at tree level (including an extra gauge link at infinity in the case of 
physical gauge [16], see also [H]), going beyond tree level requires treating lightcone singulari- 
ties [HJ[TH], associated with the x — >■ 1 endpoint. The singularity structure at x 1 is different 
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than for ordinary (integrated) distributions, giving divergences even in dimensional rcgulariza- 
tion with an infrared cut-off |19| . The singularities can be understood in terms of gauge- invariant 
cikonal-line matrix elements [19] and related to cusp anomalous dimensions (20] [2TJ [22] . 

This can be analyzed explicitly at one loop. Expansion in powers of y 2 of the coordinate- 
space matrix element ([T]) at this order gives [19l [23] 



where d is the number of space-time dimensions, /z is the dimensional-regularization scale and 
p is the infrared mass regulator. The first term in the right hand side of Eq. ([3]) corresponds 
to the case of ordinary pdfs. The lightconc singularity v 1, corresponding to the exclusive 
boundary x = 1, cancels in this term, but it is present, even at d ^ 4 and finite p, in subsequent 
terms. 

These endpoint singularities come from gluon emission at large rapidity. They imply that, 
using the matrix element ([1]), in momentum space the 1/(1 — x) factors from real emission 
probabilities do not in general combine with virtual corrections to give 1/(1 — x)+ distributions, 
but leave uncancelled divergences at fixed kj_. It is only after supplying the above matrix 
element with a regularization prescription that the distribution is well-defined. 

A possible regularization method for the endpoint is by cut-off, implemented by taking 
the eikonal line n in Eq. ([2]) to be non- lightlike [IT] [24], combined with evolution equations 
in the cut-off parameter rj = (p ■ n) 2 /n 2 [T2] [20]. Then the cut-off in x at fixed k± is of 
order 1 — x > kj_/y/4^j. Monte Carlo event generators that make use of unintcgratcd pdfs also 
implement a cut-off. We consider such applications in Sec. 3. 



An alternative method is provided by the subtractive method [251 126] , in which the direction 
n in Eq. ^ is kept lightlikc but the divergences are canceled by multiplicative, gauge-invariant 
countcrtcrms given by vacuum expectation values of eikonal operators. This method leads to 
well-prescribed counterterms |19| for the transverse momentum dependent splitting probabili- 
ties, which can be viewed as generalizing the plus-distribution regularization for kj_ =/= 0. On 
one hand, this approach has been used to relate the endpoint behavior at fixed kj_ with the cusp 
anomalous dimension |22j and investigate the role of the Mandclstam-Lcibbrandt prescription 
in lightcone gauges. On the other hand, it can likely be more useful than the cut-off to untangle 
issues of factorization and non-universality (2] [4] [27] [28] which arise beyond leading order, and 
investigate the relationship between evolution equations in rapidity and in virtuality [12 ] [20 ] [29] . 




(3) 




Figure 2: Soft gluon exchange with spectator partons. 
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Infrared subtractions analogous to those in the method [25l [26] are also discussed in the context 
of the soft-collinear effective theory [30] (under the form, however, of counterterms that are not 
automatically gauge-invariant) in the case of the Sudakov form factor [31] and of initial-state 
beam functions [32] describing the incoming jet (see also |33j). 

As noted earlier, a full treatment of factorization at the level of TMD pdfs is yet to be 
achieved [21 H]- In the hadroproduction of nearly back-to-back hadrons, factorization is bro- 
ken [4] by soft gluons exchanged between subgraphs in different collinear directions (Fig. [2]). 
(See also the analyses (34] [35] for the Drell-Yan case). The underlying dynamics is that of 
non-abelian Coulomb phase [35] , involving interactions with spectator partons [37] (which were 
treated long ago in p] limited to the case of fully inclusive Drell-Yan). It was noted in [8] that 
the factorization-breaking contributions J27j[4] are Coulomb/radiative mixing terms related to 
the contributions |38| responsible for the appearance of super-leading logarithms in di-jet cross 
sections with a gap in rapidity. 

The issue of factorization depends on developing a systematic treatment, capable of handling 
overlapping divergences in infrared regions for complex observables that involve color charges 
in both initial and final states. This motivates the subtraction techniques quoted above, in the 
version [HI [26] or in the SCET version (30j EU [32] E3J EU [35] . 

Factorization involves in general, besides TMD pdfs, a nonperturbative soft factor, also 
characterized in terms of operator matrix elements, see |25j . The first paper in Ref. |26] analyzed 
the possibility of reabsorbing the soft factor into a redefinition of the TMD pdfs (and analogous 
final-state fragmentation functions), using an explicit construction at one loop in the case of 
initial-state and final-state parton showers in DIS. More recently, an explicit procedure to 
reabsorb the soft factor into redefined TMD pdfs has been presented in Ref. [2] in the context 
of Drell-Yan. An analogous redefinition is implied by the Drell-Yan analyses in [311 ES]. It 
remains to be seen how generally this can be done. This will influence future programs of 
phcnomenological determinations of TMD pdfs, see j3J. 

2.2 TMD formulation at small x 

The case of back-to-back di-hadron or di-jet hadroproduction [H[27] illustrates that, due to the 
difficulty in disentangling soft and collinear gluon correlations between initial and final states, 
a general TMD factorization formula is still lacking. In the case of small x, however, a TMD 
factorization result holds [6] owing to the dominance of single gluon helicity at high energy. In 
this case, a TMD gluon distribution can be defined gauge-invariantly from the high-energy pole 
in physical cross sections. See EE [40j [41] [42] for recent discussions of unintegrated pdfs based 
on small x. 

The main reason why such a definition for TMD pdfs can be constructed in the high-energy 
limit is that one can relate directly (up to pcrturbative corrections) the cross section for a 
physical process, e.g. photoproduction of a heavy-quark pair [53], to an unintegrated, transverse 
momentum dependent gluon distribution (Fig. [3J). This is quite like what one does for deeply 
inelastic scattering, in the conventional parton picture, in terms of ordinary (integrated) parton 
distributions. On the other hand, the difficulties in defining a TMD distribution in the general 
case, over the whole phase space, can be associated with the fact that it is not obvious how to 
determine one such relationship for general kinematics. 

The evolution equations obeyed by TMD distributions defined from the small- a; limit are 
evolution equations in energy [5], with corrections down by logarithms of energy rather than 
powers of momentum transfer. Using this evolution, the high-energy factorization [6] can be 
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Figure 3: TMD factorization at high energy for heavy quark photoproduction. 



related, order- by-order in a s [44] , with the renormalization-group factorization [T] HH [45] • This 
allows one to describe in this framework the ultraviolet region of arbitrarily high kj_, and in 
particular rc-obtain the structure of QCD logarithmic scaling violations, see e.g. [71 14*5]. 

The above observation justifies the use of this approach for hard production physics. In 
particular, it is the basis for using Monte Carlo implementations of transverse momentum 
dependent parton showers (see e.g. [El [47] and references therein) to treat multi-scale hard 
processes at the LHC. Parton-shower applications are discussed in Sec. 3. Such applications 
have so far focused on TMD gluon distributions. Ongoing work on the generalizations needed 
to include quark channels is discussed in Sec. 4. 

Let us finally recall that extensions of the factorization results above are required if one is to 
take into account the nonlinear effects that are expected to arise in the small x region from high 
parton densities. Recent work in this direction, focusing on multiple scattering effects in dense 
targets and nuclei, may be found in [42] [48j [49] . In this respect, we note that techniques such 
as those in [5U] have been proposed to incorporate the treatment of multiple-gluon rcscattering 
graphs at small x starting from the operator matrix elements [11] [15] for parton distributions. 
They may thus be helpful for extensions to the high density region that arc aimed at retaining 
accuracy also in the treatment of contributions from high px processes. 

3 Parton-branching applications 

In this section we move to applications of the TMD formalism to parton branching methods. 
The main role of TMD splitting functions and distributions in this context is that they serve to 
take into account coherence effects of multiple gluon radiation for small longitudinal momentum 
fractions x in the initial state parton cascade. In Subsec. 3.1 we briefly recall the motivation for 
angular-ordered parton showers, and in Subsec. 3.2 we discuss the treatment of gluon coherence 
at small x and corrections to angular ordering. 

3.1 Collinear showering and soft gluon coherence 

Branching algorithms in standard shower Monte Carlo generators |51] are based on collinear evo- 
lution of the jets, both time-like and space-like, developing from the hard event. The branching 
probability is given in terms of splitting functions P and form factors A (Fig. 2]) as 
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The theoretical basis for the branching approach is the factorizability of universal splitting 
functions in QCD cross sections in the collinear limit [IJ [52], which justifies the probabilistic 
picture. 




form factors A 

Figure 4: Parton branching in terms of splitting probabilities and form factors. 

Besides small-angle, incoherent parton emission, many of the current shower generators 
also take into account further radiative contributions from emission of soft gluons, which are 
essential for realistic phenomenology |51j . To incorporate these in a probabilistic framework, one 
appeals to properties of coherence of color radiation [52] [53] [54] . Soft-gluon emission amplitudes 
factorize in terms of eikonal currents [551 156] 

n 

where pi are the emitters' momenta, q is the soft momentum, and the color charge operators Q° 
are associated with the emission of gluon a from parton i. In general, interferences are expected 
to contribute to the radiative terms relating the (n+ l)-parton process to the n-parton process. 
Nevertheless, a probabilistic branching-like picture can be recovered [57] [58] [59] by exploiting 
soft-gluon coherence. This is illustrated in Fig. [5] [47 for the case of two-gluon emission. 



a \ a 2 




Figure 5: (a) Two gluon emission from a fast quark; (b) coherence of soft gluon emission at large 
angle. 

In Fig. [5^, two soft gluons with momenta q\ and (72 are produced from a fast parton with 
momentum p. Suppose q® <C g?. We distinguish two angular regions for the softest gluon 
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qi [47]. i) When qi is at small angle from p (qi), then the amplitude can be seen as the 
sequential emission of q\ from p and of qi from p (qi). This corresponds to the standard 
brcmsstrahlung picture based on radiation cones centered around p and q\ . ii) When q2 is at 
large angle, 9 pq2 3> 9 pqi , then the directions of p and q\ can be identified and the two emission 
amplitudes act coherently to give (Fig. \Ejp) what can be seen as the sequential emission of qi 
from p and of gi from p. The reversed order of the emissions compared to case i) reflects the 
fact that the radiated gluon sees the total color charge of the emitting jet. Fig. [SJd illustrates 
that contributions of different emitters combine to give an effective contribution in which the 
emissions arc ordered in angle, so that angular ordering |57[ 1581 159) replaces energy ordering. 

The above framework of collinear showering supplemented with phase space constraints 
designed to implement the angular ordering |51j forms the basis of standard shower Monte 
Carlo event generators. We next discuss the modifications of this framework that are required 
to treat parton showers at increasingly high energies. 

3.2 Space-like parton shower at high energies 

New effects arise if one is to extend the picture of the above subsection to the case of very high 
energies, where processes with multiple hard scales become significant. The first new effect 
is that soft-gluon insertion rules [S31 [55] for n-parton scattering amplitudes can still be 

given in terms of real and virtual soft-gluon currents and 3^ v \ 

\M^ n+1 \k,p)\ 2 = |[M ( ")(fc + g,p)] t [J^] 2 M { ^(k + q,p) 

- [3^] 2 M^(k,p)} , (6) 

but the currents are modified in the high-energy, multi-scale region by terms that depend on the 
total transverse momentum transmitted down the initial-state parton decay chain [431 1601 161] . 
For this reason the physically relevant distribution to describe space- like showers at high energies 
is not an ordinary parton density but rather a TMD parton density. 

The next point concerns the structure of virtual corrections. Besides Sudakov form-factor 
contributions included in standard shower algorithms |51[ 152) . one needs in general virtual- 
graph terms to be incorporated in transverse-momentum dependent (but universal) splitting 
functions |60j . These allow one to take account of gluon coherence not only for collincar-ordercd 
emissions but also in the non-ordered region that opens up at high ^/s/p±, where y/~s is the 
total center-of-mass energy and p± is the typical transverse momentum of a produced jet. 




Figure 6: (left) Coherent radiation in the space-like parton shower for x <C 1; (right) the unintegrated 
splitting function V. including small-x virtual corrections. 
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The resulting structure of the parton branching, depicted in Fig. El differs from that in 
Eq. ((4]): the branching probability is k^-dependent, and part of the virtual corrections are 
associated to the (unintegratcd) splitting functions. Schematically one has, using the recursion 
relation ©, 

Q(x,k±,fi) = Qo(x,k±,fi) + J ^ J ^jj- &(p- zq) 

X A(^zq)V(z,q,k ± )g{-,k ± + (l-z)q,q) , (7) 

z 

where Q is the unintegrated gluon distribution, A is the form factor, and V is the unintegrated 
splitting function (Fig. [6]) . The kernels V depend on transverse momenta and include part of 
the virtual corrections, in such a way as to avoid double counting with the Sudakov form factor 
A, while reconstructing color coherence not only at large x but also at small x in the angular 
region (Fig. 

a/x > ot\ > a , (8) 

where the angles a for the partons radiated from the initial-state shower are taken with respect 
to the initial beam jet direction, and increase with increasing off-shellness. 

In terms of high-energy power-counting, the effects of region (|HJ) are potentially enhanced 
by terms 

a^ln fe+m ^/pj. ■ (9) 
In inclusive processes, coherence leads to strong cancellations between reals and virtuals so 
that terms with m > 1 in Eq. © drop out, and high-energy corrections are at most single- 
logarithmic [44, 62, 63 . For exclusive jet distributions such cancellations are not present and 
one may expect stronger enhancements. The implementation of coherent-branching effects 
associated with high-energy logarithms is, from the point of view of jet physics, the main 
motivation for developing the formalism of unintegrated parton distributions and implementing 
it in shower Monte Carlo event generators. See [M] for recent phenomenological investigations 
of such effects on the structure of angular correlations and multiplicity distributions in multi-jet 
final states, and |65) for applications to the jet structure associated with heavy mass states. 

Monte Carlo implementations of the modified branching described in this section include 
the parton-shower event generators Cascade [66], Ldcmc [67]. See [H[47l[68] for further refer- 
ences. These implementations take into account gluon coherence effects and include TMD gluon 
distributions. An extension to valence quark distributions is given in [69) . In the next section 
we describe recent work to go beyond this approximation and include sea quark contributions. 



4 Sea quark distribution and Drell Yan production 

Monte Carlo calculations based on TMD approaches have so far been carried out mostly in a 
quenched approximation in which only gluon and valence quarks arc included. In terms of the 
parton branching kernels, this corresponds to taking into account only contributions from the 
splitting vertices (a) and (b) in Fig. [7J In this section we briefly describe the work [§] to go 
beyond this quenched approximation and treat TMD sea quark contributions. The main focus 
of this work is to take into account effects of the splitting process in Fig. (7{c) at TMD level. 

Let us recall that early attempts [70J [7T] to treat the unintegrated pdf evolution beyond 
the quenched approximation include quarks via splitting probabilities to lowest order of per- 
turbation theory, neglecting any transverse momentum dependence in the branching. In |72j 
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k^-dcpendent kinematic corrections are included, while the splitting kernels are still taken in 
lowest order. 




(a) (b) (c) (d) 



TMD gluon and 
valence quark 



Quark emission 
contribution to TMD sea 



Figure 7: Splitting vertices for TMD pdf evolution. 

Also, a program to perform shower Monte Carlo evolution at unintegrated level has recently 
been proposed [73] based on the expansion [45] in two-particle irreducible (2PI) kernels. This 
program is formulated at the next-to-leading order. However, this is limited at present to the 
flavor nonsinglet sector, i.e. does not include small- a; logarithmic corrections which are present 
beyond NLO in flavor singlet distributions. 

The approach taken in [3] is based on the high-energy form of the (off-shell) TMD quark 
Green function defined in |44j . This is obtained by generalizing to finite transverse momenta, 
in the high-energy region, the 2PI expansion [45] . The k^-dependent gluon-to-quark splitting 
kernel is depicted in Fig. [SJb) . Thus, unlike the approaches in [701 [13 [Z2, on one hand, 
and in [73], on the other hand, Ref. [9] constructs TMD quark distributions by using the kj_- 
dependent gluon-to-quark splitting function, which governs sea quark evolution to all orders 
in perturbation theory in the small- a; limit. Refs. [§] I74| investigate the phcnomcnological 
relevance of the finite k± terms for Drcll-Yan production (Fig. [Sfa)). 

Schematically, the TMD gluon-to-quark kernel has the form 

V g -* q {z;qx,k ± )=P£K*) (l + £ K(z)(kl/ q i)A , (10) 

where is the lowest-order DGLAP splitting function, z is the longitudinal momentum 
transfer, fcj_ and q± are respectively the gluon and quark transverse momenta in Fig.[5](b), and 
all coefficients b n (z) of the finite terms in k± are known from [33]. Although it is evaluated 
off-shell, the splitting probability in Eq. ([TO]) is universal [44] [75]. The approaches in [70] 
and [7TJE2] only include the term P'°) in Eq. (p~0|) . while Ref. [5] includes the full series. The 
finite terms in k± are the ones responsible for determining small- a; logarithmic corrections to 
flavor-singlet quark evolution to all orders in a s . 

In Ref. [TUJ contributions from the TMD flavor-singlet quark distribution are implemented 
in the parton shower Monte Carlo event generator Cascade [5B]. This constitutes a starting 
point to systematically include quark emissions in the parton shower at TMD level. At present 
the implementation is done in such a way that the shower couples to quarks only once, while a 
complete implementation will allow this to occur arbitrarily many times. Nevertheless, already 
the present implementation incorporates for the first time the small- a; dynamics encoded in 
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(b) 

Figure 8: (a) qq Drell-Yan production; (b) g — > q splitting contribution to sea quark distribution. 

Eq. (|10[) in a parton shower, and makes it possible to treat in this framework hard processes 
induced by sea quarks on the same footing as processes induced by gluons. 

Refs. [IH] present results within this framework for Drell Yan production processes. This 
is phcnomenologically relevant, because Drell Yan production at the LHC (Fig. \E)(a)) receives 
large contributions from sea quark scattering at small x [68j . Such contributions affect many 
aspects of LHC physics, as Drell-Yan processes are instrumental in precision electroweak mea- 
surements, in luminosity monitoring and pdf determinations, and in new physics searches. 
Results on high-energy Drell-Yan at TMD level have so far been obtained in the qg* chan- 
nel [7B|, and in the associated production channel Z/W+ heavy quarks [77]. The qg* channel 
is of direct phenomenological relevance for forward Drell-Yan production. Ref. [9] evaluates the 
perturbative coefficients for the coupling of the TMD sea quark distribution to Drell-Yan by 
using the "reggeized quark" calculus [75] . This extends the high-energy effective action formal- 
ism [7S] , currently explored at NLO [HI] , to amplitudes with quark ^-channel exchange [ST] [52] . 
In [S] [TU] the method is employed to investigate predictions for Drell-Yan production in the 
forward region (Fig. [9]). 




Figure 9: Forward Drell-Yan production. 
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In Figs. 1101111 we report numerical results for Z-boson rapidity and transverse momentum 
spectra. We present results for two possible choices of the scale fi in the shower evolution 
equation (J7J: one in which [i is set equal to the hard scale, defined by 

M 2 =pi + Ml , (11) 

where p± and Mz are the transverse momentum and mass of the Z boson; another in which [i 
is set equal to the maximum-angle scale determined by the angular ordering kinematics, given 
by 

, _ ql + {l-z)kl 
A* - (1-3)2 > W 

where the variables are as specified below Eq. (jlOp . Since we are working in the forward region, 
we also include results obtained from the qg* matrix element. 



transverse momentum Z-boson 
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Figure 10: Z-boson transverse-momentum and rapidity spectra. 

All three curves in Figs. 1101111 correspond to calculations which are applicable at forward 
rapidity. The shape of the rapidity spectrum in Fig. 1 1 01 falling off for central production, can 
thus be regarded as defining the kinematic region of applicability. The transverse momentum 
spectrum in Fig. 1101 indicates that the three calculations differ at low p t , but they all converge 
for large enough p t . This is illustrated in more detail in Fig. 1111 The difference between the ma- 
trix element and factorized calculations in the small-pt region is mostly due to quark s-channel 
contributions. This region is dominated by Sudakov form factor effects, giving the turn-over in 
the p t spectrum, which are sensitive to the different choices of scales in Eqs. (|ll|) . (fT!?l) . The dif- 
ference between the curves disappears in the large p t region as the quark t-channel contribution 
dominates, and contributions from the parton showers also become relatively insensitive to the 
different scale choices, because transverse momentum ordering sets in driven by the high p t . 

The numerical results in Figs. 1101111 can be seen as a consistency check on the approach 
based on TMD parton branching as regards forward Drcll-Yan production. The result [9] for the 
unintegrated flavor-singlet quark distribution is however more general, and is a necessary ingre- 
dient to develop a description of vector boson production and associated final states throughout 
the rapidity range of forward and central regions. This is the subject of forthcoming work. 
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Figure 11: Numerical comparison in the small-pT and large-pr regions. 



5 Summary 

We have discussed recent progress on transverse momentum dependent parton distributions. 
Relying on TMD factorization results at small x, we have discussed applications to branching 
Monte Carlo methods for initial-state parton showers in high-energy collisions. Calculations 
based on these methods have mostly been developed within a quenched approximation in which 
only gluon and valence quark effects are taken into account at unintegrated level. We have 
reported work to go beyond this approximation by including sea quark contributions. 

The method is based on taking into account flavor-singlet quark evolution via the transverse 
momentum dependent gluon-to-quark splitting kernel in Eq. (|T0|) . We have used this along 
with the Monte Carlo implementation [66] of the small- a; coherent branching equation given in 
Eq. ([7]). Compared to previous approaches to treat the TMD sea quark distribution, the main 
feature of this approach is that it includes, in addition to the lowest-order splitting function, the 
full series of finite-fc^ terms in the gluon-to-quark two-particle irreducible kernel. This allows 
one to sum small- cc logarithmic corrections to flavor-singlet observables to all orders in a s . 

We have used this framework to obtain numerical predictions for forward Z-boson produc- 
tion. To this end, the pcrturbative coefficients for coupling the TMD sea quark distribution to 
Drcll-Yan production have been determined by using the formalism |78j for high-energy quark 
i-channel exchange. We have investigated the dependence of the results on the shower evolution 
scale. 

The method proposed in this work is implemented at present by letting quarks interact 
with the shower only once. Using the ingredients discussed in this work, this could however 
be extended by constructing a fully coupled shower evolution in the flavor-singlet sector. Also, 
only one coupling of the sea quark to vector bosons is evaluated at TMD level at present. 
This should be extended in order to treat Drcll-Yan throughout the rapidity range of forward 
and central regions. Nevertheless, the results presented here, being the first implementation of 
small- a; sea-quark dynamics in a parton shower, constitute a starting point to include quark 
emissions systematically and to treat hard processes initiated by sea quarks on the same footing 
as processes initiated by gluons. 
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